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Abstract 

a 

A mesh condition is developed for linear finite clement approximations of anisotropic diffusion- 
convcction-rcaction problems to satisfy a discrete maximum principle. Loosely speaking, the 
condition requires that the mesh be simplicial and O(||6||oo^ + ||c|| 00 /i 2 )-nonobtuse when the 
dihedral angles are measured in the metric specified by the inverse of the diffusion matrix, 
where h denotes the mesh size and b and c are the coefficients of the convection and reac- 
tion terms. In two dimensions, the condition can be replaced by a weaker mesh condition (an 
0(||b||cx>/i+ ||c||oo/i 2 ) perturbation of a generalized Delaunay condition). These results include 
many existing mesh conditions as special cases. Numerical results are presented to verify the 
theoretical findings. 
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^ 1 Introduction 



O 

(N 



We are concerned with the linear finite element (FEM) solution of the anisotropic diffusion equation 

-V ■ (OVu) + b-Vu + cu = f, in n (1) 



subject to the Dirichlet boundary condition 

H 



u = g, on 8ft (2) 

where Q C R d (d > 1) is a connected polyhedron and D = B(a?) 6 M. dxd (the diffusion matrix), 
b = b(x) € M d , c = c(x), f = f(x), and g = g(x) are given, sufficiently smooth functions defined 
on Q. We assume that for any H>(x) is symmetric and strictly positive definite and functions 

b and c satisfy 

c(x) - • b(x) > 0, c(x) > 0, Vxe n. (3) 
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It is known (e.g., see [8]) that the solution of boundary value problem (BVP) ([I]) and ([2]) satisfies 
the maximum principle. 

The numerical solution of BVP ([I]) and ([2]) has attracted considerable attention from scientists 
and engineers. The BVP is a prototype model for anisotropic diffusion problems which arise in 
various fields such as plasma physics [TUl [Til ES] ; petroleum reservoir simulation [U [23] , and image 
processing |24[ 129j. Moreover, it has been amply demonstrated that a standard numerical method, 
such as a finite element, a finite difference, or a finite volume method, does not necessarily satisfy 
a discrete maximum principle (DMP) and may produce unphysical solutions that typically contain 
spurious oscillations, undershoots, and overshoots. Furthermore, designing a numerical scheme to 
preserve the maximum principle is an important research topic in its own right. As a matter of fact, 
considerable work has been done in the past to develop numerical schemes to satisfy DMP; e.g., 
see [21 El HI El EH UH HS1 HSl HE EH EH [30] for isotropic diffusion problems (D = a{x)I with a(x) 
being a scalar function) and |Zl HD1 HI1 HH HS1 H3 QS1 HD1 12H 1221 US] for anisotropic diffusion 
problems. In particular, it is shown in [6] that the linear FEM satisfies DMP when the mesh is 
simplicial and satisfies the so-called non-obtuse angle condition which requires that the dihedral 
angles of all mesh elements be non-obtuse. In two dimensions the condition can be replaced by 
a weaker condition (the Delaunay condition) which requires that the sum of any pair of angles 
opposite a common edge is less than or equal to tt |27j . Similar results have been obtained recently 
for anisotropic diffusion problems in |12[ [19] . 

It is pointed out that most of the existing work has been concerned with problems without 
convection terms. For continuous problems, it is known (e.g., see jS]) that convection terms have 
no effect on the satisfaction of the maximum principle by the solution. However, the situation is 
different for discrete schemes. The main difficulty comes from the fact that discrete convection terms 
typically do not vanish at an interior maximum point and the entries of the corresponding matrix 
can be both positive and negative. A few researchers have tried to address the issue for isotropic 
diffusion problems. For example, Xu and Zikatanov [30] employ a special number treatment for 
convection terms so that they have no effect on the DMP satisfaction by the discrete solution. 
Burman and Ern [3J propose a nonlinear stabilized Galerkin approximation of the Laplace operator 
which satisfies DMP on arbitrary meshes and for arbitrary space dimension without resorting to the 
non-obtuse angle condition. They prove that the result can extend to diffusion-convection-reaction 
problems with constant diffusion coefficient when the mesh is locally quasi-uniform. More recently, 
Wang and Zhang [28j study quasilinear isotropic diffusion-convection-reaction problems and show 
that linear finite element approximations satisfy DMP when the mesh is 0(||6||oo^+ ||c|| 00 /i 2 )-acute 
(i.e., the dihedral angles of all mesh elements are less than or equal to | — 711161100/1 — 72HCHOO/1 2 
for some positive numbers 71 and 72). On the other hand, no work has been done for anisotropic 
diffusion-convection-reaction problems. 

The objective of this paper is to develop a mesh condition for linear finite element approxima- 
tions of anisotropic diffusion-convection-reaction problems ([!]) and (|2j) in any dimension to satisfy 
a discrete maximum principle. We shall use the approach of |19] to show the stiffness matrix 
associated with the linear finite element discretization to be an M-matrix and have non-negative 
row sums, with the focus on the treatments of the convection and reaction terms. We shall also 
investigate the two dimensional case where a weaker sufficient condition can be developed. 

The paper is organized as follows. A linear finite element discretization for BVP (TlJ) and Q is 
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introduced in Section 2. In Section 3 geometric properties of the gradient of linear basis functions 
are studied. A general mesh condition valid in any dimension and a specific and weaker condition 
in two dimensions are developed in Section 4, followed by numerical results in Section 5. Finally, 
Section 6 contains conclusions. 



2 Linear finite element formulation 

We consider the linear finite element solution of BVP ([!]) and ^ . Assume that an affine family of 
simplicial meshes {Th} is given for f2. Let 

U g = {v£H 1 (n) \v\ an =g}. 

Denote by U^ h the linear finite element space associated with mesh Th, where g h is a piecewise 
linear approximation to g on the boundary. A linear finite element approximation u h 6 U^ h to 
BVP and (g]) is defined by 

/ (V/) r BV/di+ / v h (b ■ Vu h )dx + [ cu h v h dx= [ fv h dx, \/v h eU{}. (4) 
Jn J n Jn Jn 

The above equation can be rewritten as 

\K\ (Vv h fO K Vu h + [ v h (b-Vu h )dx 

KeT h KeT h Jk 

+ J2 [ cuhvhdx = Yl [ f yhdx i Vv h eUi} (5) 
KeT h K Ken K 

where \K\ is the volume of element K and Dk is the integral average of D over K, viz., 



1 

W\ 



1 dx. (6) 

K 



Scheme ^ can be expressed in a matrix form. Denote the numbers of the elements, vertices, 
and interior vertices of mesh Th by N, N v , and N v i, respectively. Assume that the vertices are 
ordered in such a way that the first N V i vertices are the interior vertices. Then Uq and u h can be 
expressed as 

Uq = span{^i, • • • , 4> Nvi }, (7) 
N vi N v 

j=l j=N vi +l 

where <f>j denotes the linear basis function associated with the j-th vertex, aj. The boundary 
condition ([2]) is approximated by 



u 



.j=g{aj), j = N vi + 1,...,N V . (9) 



Substituting ^ into ([5|, taking v h = 4>j [j = l,...,N v i), and combining the resulting equations 
with ([9]), we obtain the linear algebraic system 

Au = f, (10) 
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where u = (u\, ...,u Nvi ,u Nvi+ i, ...,u Nv ) T , f = (/i, fN vi , 9N vi +l, 9n v ) T , 



A 



An An 
/ 



(11) 



and I is the identity matrix of size (N v — N V i). The entries of the stiffness matrix A and the 
right-hand-side vector / are given by 

Oij= Y \K\(V4>i) T B K V^3+ Y / 4>i{b-V(j) j )dx+ Y / cfafodx, (12) 
KeT h Kan K Kan K 

i = 1, ...,N vi , j = 1, ...,N V 



fi= Y f<l>i dx i i = l,-,N vi . 



(13) 



Ken 



In the following sections we shall investigate under what condition on the mesh the solution of ^ 
satisfies a maximum principle. A key to this investigation is to understand geometric properties of 
the gradient of linear basis functions which are to be described in the next section. 



3 Geometric properties of the gradient of linear basis functions 

Let K be an arbitrary simplex with vertices oi, a<i, o-d+i- Denote the face opposite to vertex aj 
(i.e. the face not having aj as its vertex) by Si and its unit inward (pointing to aj) normal by nj. 
The distance (or height) from vertex aj to face Si is denoted by hi. The result of the following 
lemma exists in literature; e.g., see [2"1 [T^l [50]. 

Lemma 3.1 For any simplex K € M. d , the gradient of linear basis function (pi associated any 
vertex ai (i = 1, d + 1) is given by 

VJ>i = £. (14) 



It is remarked that Brandts et al. [2] have obtained the same result using the so-called q- vectors 
defined through the edge matrix of elements. Specifically, they show that Qj, a q- vector associated 
with face Si, is an inward normal to Si, has the length 1/hi, and is equal to V0j; i.e., 

Qi = T-n, = V0j, i = l,...,d+l. (15) 
i%i 

These q-vectors will be used frequently in the remaining of the paper. 

The next property of gradient of linear basis functions is related to the diffusion term in stiffness 
matrix (12) for the case B>k = I. Denote the dihedral angle between any two faces Si and Sj (i 7^ j) 
by ocij. It can be calculated as the supplement of the angle between the inward normals to the 
faces, i.e., 

cos(ajj) = - m ■ rij = - - — 1 1 , i / j. (16) 

9< h II II 



(In fact, (16) is often used as the definition of the dihedral angle.) A sketch of the q-vectors, 
dihedral angles, and heights of an element are shown in Fig. [TJ 

The result of the following lemma is also known in literature; e.g., see [2"1 l9l [12]. 
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Figure 1: A sketch of unit inward normals, dihedral angles, and heights of element K. 



Lemma 3.2 For any simplex K £ M. d , we have 

\K\ 



It reduces to 



in two dimensions. 



\K\(V<t>iYV^ 



hihj 



cos(ay), i^j. 



cot(ajj), i^j 



(17) 
(18) 



Proof. Equation (17) follows from Lemma 3.1 and (16) 



In two dimensions, K is a triangle. Consider the case with % = 1 and j = 2. From Fig. [TJ we 
have 

hi hih 2 



\ K \ = y H a 2 - a 3l 



2sin(ai 2 ) ' 



Combining this result and (17) gives (|18j). □ 
We now study the diffusion term \K\(\7(pi) T Il>K^4 > j f° r general symmetric and positive definite 



matrix D# using Lemma 3.2 Define 



G K {x) 



x : K -> K, 



(19) 



where K = G(K). Obviously, K is also a simplex in M. d . For any vertex otj, we denote the 
corresponding vertex, face, height, and q-vector of K by Oj, Si, hi, and respectively. We have 



a, 



=det(B K )-2\K\, q t 



(20) 



where || • ||d k denotes the distance measured in the metric Hk- The derivations of the first three 
relations are trivial. To derive the last two, we first notice that 

i 

4>i(x) = cj)i(p 2 K x) = (j>i(x). 



Then from (15) we have 



V(j)i 



which gives the second last relation in (20). The last relation is obtained by taking the norm of the 
above equation. 



5 



To obtain the relation between hi and hi, we rewrite the last relation in (20) as 

hi 



from which we obtain 

hi ~ hi , 

—j=^= <K< , (21) 

V /x inax 

where A max (B/f ) and A m j n (Dx) denote the maximum and minimum eigenvalues of Bjf, respectively. 

_ _ _x _ 

Denote the dihedral angle between faces Si and Sj by ay D -i. Since Si = D K 2 Si and Sj = 

_ i 

D^- 2 Sj , it can also be viewed as the dihedral angle between Si and Sj measured in the metric D^- 1 . 



Moreover, from ( 16 ) we see that the angle can be calculated by 



cos (%',D^ ) = " n~ ii n4 H = " TTTn i^TJT- • ( 22 ) 



" K ' \\Qi\\ ■ \\Qj\\ \\Qi\\n K \\Qj \\to k 



We now go back to the quantity \K\(SJ 4>i) T I5) (f>j ■ Notice that 



Applying Lemma 3.2 to K and using relations (20), we have the following lemma 



Lemma 3.3 For any simplex JfeR and any symmetric and positive definite matrix Ok, we 
have 



{KKV^DkV^ = - Jg! yj *) 2 cos(a B -0, i^j. (23) 



It reduces to 



in two dimensions. 



|ir|(V<&) T PieV^ = _ det (^) 2 cot( % . B -i), »Vj (24) 



4 Mesh conditions for DMP satisfaction 

In this section we study the mesh conditions under which the linear finite element scheme ^ 



satisfies DMP. The main conclusions are given in Theorems 4.1 and 4.2 
Theorem 4.1 If the mesh satisfies 

h^f ll&lloo k hfhf II ell oo k . f , 

+ ; hi s - 73 \7f, o\ <cos(a t , B -i), (25) 



Xmin(0 K ) (d+1) X min (D K ) (d + l)(d + 2)- V W 



where ||b||oo,A: = rnaxa; e # ||b(a;)||, HcHoo^ = m&x X £Kc(x), and hf 's and a i<jn -i 's are the heights 

' ^ F— 1 

and dihedral angles of element K, respectively, then the linear finite element scheme ([5p for BVP 
pi) and pi) sate/ies £>MP. 
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Proof. Following [19] we prove this theorem by showing that stiffness matrix A defined in 
(11) and (12) has non-negative row sums and is an M-matrhPJ From Stoyan |26|, Theorem 1], this 
implies that scheme ([5]) satisfies DMP. 

(1) We first show that matrix A has non-negative row sums. Notice that we only need to show 
the first N v i row sums are non-negative. Using the fact Ylf=i $j{ x ) = 1 an ^ the assumption c > 
(cf. Q), from (12) we have, for i = 1, ...,N v i, 

N v I N v \ / / N v \ \ 

^ aij = \ R \ ( V ^) T Bk v + E / 0, kv &] )dx 

i=i A'eTh \j=i J KeT h Jk \ \j=i 



+ E I c<f>i \J2<l>j) dx 
Ken Jk \i=i / 

2 c 4>idx 



KeT h 

> 0. (26) 

(2) Next, we show that A is a Z-matrix; i.e., 

at, < 0, Vi^j,i = l,...,N vi , j = l,...,N v (27) 
an > 0, i = l,...,N vi . (28) 

Recall from Ciarlet [3 Page 201] that 

li^i r , , . lifi 



M * = dn' y ^" = (d + i)(^ +2 ) ' (29) 

Keuji KeuJiHoOj 

where cjj and are the element patches associated with vertices Gtj and aj, respectively. We have 

a,ij= V ( \K\ (X7<fti) T H>k V0j + / <j>i(b-V(t>j)dx+ / c fa <j)jdx) (from (Q) 
KeuiC\Wj V ^ ' ^ / 

< (V4>i) T B> K V(t)j + 4>i\b-nf\dx+ / c ^ 0j da; ) (Lemma pDj 
Ke*^ V ^ ^ ^ / 

< E (l*l CV*f V* + ™^ + T^^y) (from @ ) 

= > — ~ — ~ — COS a.. n -iH it-, ; — r t~, n~, r Lemma 6.6) 

K^A h?h? h*(d + l) (d+l)(d + 2)) LJ 

V" \ K \ (_ h?hf\Moo,K hfhf\\c\\^ K \ 

K\ ( , , hf llblU.jf /if /if llclU^ \ 



I v y X min {B K )(d + l) X min (B K )(d+l)(d + 2) r 



(from (21)) 



1 Matrix yl is called an M-matrix if it is a Z-matrix (see (27 1 and (28 1 below) and satisfies A 1 > (i.e., all entries 
of its inverse are nonnegative) . 
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Combining this with (25) implies (27). 



On the other hand, for i = 1, ...,N v i, 

a u = ^ \K\ {V^f B K V^i+ I <f>i(b-V<f>i)dx+ ! ctfdx (from @) 

> / 4>i{b-V(j)i)dx+ / c^dcc 
Jf7 Jn 

= I (c V • b)<p 2 dx. (Gauss' divergence thm) 

The assumption ^ implies that a« > 0. Thus, stiffness matrix j4 is a Z-matrix. 

(3) We now show that An, the northwest block of matrix A, is an M- matrix. This is done by 
showing An is positive definite. For any vector v = (vi,V2, ...,vn vi ) T , we define v h = v i^i ^ 
Uq. Notice that Vv h is constant on K. As in the proof for an > 0, from (12) we have 



v T A u v = Yl \K\ (Vv h f B K Vv h + [ v h (b-Vv h )dx + I c(v h ) 2 dx 
> \K\ {Vv h ) T B K Vv h + ! (c - • b){v h ) 2 dx > 0. 



KeT h 

Moreover, from the above inequality, v T Anv = implies v = constant, which in turn implies 
v h = due to the fact that v h G Uq. From these, we know that An is positive definite. Since An 
is a Z-matrix, so it is an M-matrix. 

(4) Finally, we show matrix A is an M-matrix by showing the inverse of A is positive. From 



( 11 ), the inverse of A is given by 



A' 



A n * -A^Au 
/ 



(30) 



Using the fact that Aq 1 > and A12 < 0, then A 1 > 0, which, together with the fact that A is a 
Z-matrix, implies that A is an M-matrix. □ 

It is remarked that the result in the above theorem is consistent with those obtained by Brandts 
et al. [2] for (anisotropic) diffusion-reaction problems and Wang and Zhang [28J for (anisotropic) 



diffusion-convection-reaction problems. Loosely speaking, (25) requires 



or 



cos(a.. D -i) > 0(h\\b\\oo) + O^lcHoo) 



< C^-jj-l < I - 0(/l||6||oo) - 0(h 2 |lc]|oo) 



(31) 



(32) 



for all dihedral angles, where h = max hu is the maximum element size. In other words, if the 

KeT h 

mesh is 0(h)-acute in the metric D _1 for the case b ^ or 0(h 2 )-acute in the metric D^ 1 for the 
case b = and c ^ 0, then the linear finite element solution of |7p and |lp satisfies a DMP. 

It is known that the acute or nonobtuse angle condition can be replaced by the weaker, so-called 
Delaunay condition in two dimensions for a linear finite element solution to satisfy a DMP; e.g., see 
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Strang and Fix |27j for the anisotropic diffusion case and Huang |12j for the anisotropic diffusion 
case. In the current situation with convection and reaction terms, a similar weaker condition can 



also be obtained in two dimensions. The argument is almost the same as that of Theorem 4.1 
except that Step (2) of the proof needs to be fine-tuned. Let eij be the edge connecting vertices 
aj and aj (i = 1, ...,N v i, j = 1, ...,N, i ^ j). Denote the two elements sharing ey by K and K' . 
From Step (2), we have 

an < \K\(y^\ K fB K V^\ K + lK l m °°> K + |iV1 IHk ' K 



hf(d+l) (d+ l)(d + 2) 



+ \K'\{V(t)i\K') T ^K^HK' + 



I if' I \\b\\oo,K' . \K'\ Hclloo^/ 



hf(d+l) (d + l)(d + 2) 



det(D^)a 



cot(a- 



det(D jftr /)2 



cot (a 



\K\ ||b||oo,x \K\ llclloo^ ^ \K'\ Hblloo,^ |if'| llclloo,^' 



fcf(d + l) (d + l)(d + 2) fcf(d+l) (d+l)(d + 2)' 



(Lemma 3.3 ) 
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where ct-- n -i and a- n-i are the angles of K and K' respectively, that face the common edge en. 
From this we can conclude that the linear finite element solution in 2D satisfies a DMP if the mesh 
satisfies 

I if I ||&||oo,a: \K\ HcIIoo^ \K'\ ||6|| 00) ji'/ \K'\ Wc^^ 1 
hf(d+l) (d+l){d + 2) + hf'(d + l) + (d+l)(d + 2) 



< 



det(I 



cot (a- 



■ | det(D^)2 
-i ) + cot(a WJl -i ( 



for all internal edges. Following [12], we can rewrite the above inequality as 



(33) 



where 



< 



1 



"ij.O- 1 + a ij,D-) + arCCot 



det(D^) 2C(K,K',j) 
det(B K ) COl ^' K V J de t(D K ) 



+ arccot 



det( 



"ATJ 



C(K,K',j) 



y det(DA-/) 
|if| ||6||oo,K 



cot (a. 



2C(K,K',j) 
• D ^ lj V det(D^) 



< 7T, 



, |if| ||c||oo,X , \K'\ ||b||oo,A-' I -K"' I llcllocJC' 



/if(d+l) (d+l)(d + 2) (d + l)(d + 2)' 



The following theorem summarizes the above analysis. 



(34) 



(35) 



Theorem 4.2 If (34) holds for all internal edges of the simplicial mesh 7~h, then the linear 
finite element scheme |5|) for BVP and Q) in two dimensions satisfies DMP. 



Loosely speaking, (|34|) can be written as 
1 



< 



^B- 1 + + arCCOt 



det(D X ' 



arccot 



det(DA-) 

— -— r COtia- ■ n -l 



det(D^ 

<TT-0(h\\b 



0(h 2 \\c\\ 



(36) 
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For the case where D = I, b = 0, and c = 0, it is easy to see that (36) reduces to the Delaunay 
condition: a%j k + o-ijK' < tt- Moreover, for the case without convection and reaction terms, (36) 
gives the Delaunay-type mesh condition obtained by Huang for two dimensional anisotropic 
diffusion problems. 

5 Numerical examples 

In this section we present numerical results obtained for four 2D examples to verify the mesh 



condition (25) and (36). In all but Example 5.4 the convection vector b is taken as a constant 
vector with equal, positive x and y components, i.e., b = ^^(l, 1) T . 

Example 5.1 The first example is in the form of ([I]) and Q, and the coefficients are given as 

c=0, / = 0, g(x,0)=g(16,y) = 0, 

j 0.5y, for0<y<2 / 1, for < z < 14 

q(U,y) = < and g(x, lb = < 

v ' | 1, for2<y<16 x ' 1 8 - 0.5a;, for 14 < x < 16. 

For this example, the diffusion matrix is taken as the identity matrix, i.e., D = I. This is an 
isotropic homogeneous diffusion problem. Note that the example satisfy the maximum principle 
and their solutions stay between and 1. 

An acute-type mesh is used in the computation. Such a mesh is obtained by partitioning each 
square element of a uniform mesh into eight triangles with acute angles; see Fig. [2] The maximum 



angle of the mesh is 0.497T and thus condition (25) holds when the mesh size is sufficiently small. 



Fig. [3] shows the contours of the linear finite element solutions obtained for N = 9800 and 
N = 20000. (N is the number of elements.) There are no undershoot nor overshoot for N = 20000 
whereas both undershoots and overshoots occur for the case with N = 9800. In Fig.ga), -u min is 
shown as functions of the number of elements N . From the figure one can see that —u m i n decreases 
as the mesh is refined and the decrease rate is about quadratic initially and then exponential near 
N = 10000. Moreover, 

~ u min becomes zero (more precisely, at the level of roundoff error) after 
around N = 17000. This is consistent with Theorem 14.11 which states that there are no undershoot 
nor overshoot when the mesh size is sufficiently small. 



To further verify Theorem 4.1, we fix the number of elements at N = 3200 and let ||&||oo vary. 
Quantity —u m i n is plotted in Fig. [4](b) as a function of ||6||oo- From the figure, we can see that 
there is no undershoot until ||6||oo ~ 4. Then —u m in increases exponentially until ||6||oo ~ 20 where 
the increase rate is about linear as ||6||oo increases. 

Finally, it is pointed out that a similar behavior can be observed for the overshoot. The results 
are omitted here to save space. □ 



Example 5.2 In the second example, the BVP ([!]) and ([2]) with all the coefficients being the 



same with Example 5.1 except the diffusion matrix is used. The diffusion matrix is taken as 



(x,y) 



500.5 499.5 
499.5 500.5 



10 



2 4 6 



Figure 2: A typical mesh (N = 200) used for Example 5.1 
(a): N = 9800 (b): N = 20000 
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y ./// 
// / 
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Figure 3: Contours of the linear finite element solutions for Example |5.1 



This matrix represents a homogeneous but highly anisotropic diffusion process. A mesh with right 
triangle elements (see Fig. [5]) is used for this example. Such a mesh is obtained by dividing each 
square element of a uniform mesh into two right triangular elements. Although each element of 
the mesh is a right triangle (in the Euclidean sense), the maximum angle is 0.49-7T when measured 
in metric D _1 . Thus, the mesh is of acute-type in the metric and condition (25) can be satisfied if 
the mesh size is sufficiently small. 

Contours of linear finite element solutions are shown in Fig. [6] while the undershoot is plotted 
as functions of iV and ||b||oo i n Fig. [7J From these results we can observe a similar behavior of the 
undershoot and overshoot as in Example |5.1[ i.e., they occur only for relatively coarse meshes or 
relatively large ||b||oo- The behavior is consistent with Theorem 4.1 n 



Example 5.3 In this example, the same BVP Q and ([2]) with Example 5.1 except the diffusion 
matrix is used for this example. The diffusion matrix is taken as 



50 12 
12 50 
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(a): b= (10, lO) 5 



(b): N = 3200 





100 1000 10000 100000 

N: numbers of elements 



Figure 4: The undershoot, —u m i n , is plotted as a function of the number of elements iV in (a) and 



as a function of ||6||oo in (b) for Example 5.1 




Figure 5: A typical mesh (N = 200) used for Example 5.2 (with anisotropic 



This diffusion matrix has a weaker anisotropy than that in the previous example. 

The mesh used for Example |5.1| (see Fig. [2J is also used for this example. Recall that the mesh 
is acute in the Euclidean sense. When measured in metric D , however, the maximum angle of 
the mesh is 0.557T and the maximum sum of any pair of angles opposite a common edge is 0.977T. 



Thus, the mesh will satisfy (36) but not (25) when its size is sufficiently small. 



Contours of numerical solutions are shown in Fig. [8] while the undershoot is plotted as functions 
of N and ||&||oo in Fig. [9j A similar behavior of the undershoot and overshoot can be observed as 
for the two previous examples. n 



Example 5.4 In the last example, we consider the BVP (111 and S on domain = [0, l] 2 \[g, |] 2 



with 



1000(y- 0.5), 1000(a: -0.5) , c = 0, / = 0, g = on dQ out , g = 2 on d£l in , 



where dQ on t and <9Sl; n are the outer and inner boundaries of 0, respectively. The diffusion matrix 
is taken as 

cos a — sin a \ I 1000 W cos a sin a \ 
sin a cos a \ 1/1— sin a cos a I 
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(a): N = 450 



(b): N = 7200 





Figure 6: Contours of linear finite element solutions for Example 5.2 
(a): b = [200, 200] T (b): N = 3200 
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Figure 7: The undershoot, —u m i n , is plotted as a function of the number of elements TV in (a) and 
as a function of 



in (b) for Example 5.2 



where a = ir sin(x) cos(y). This diffusion matrix is anisotropic and heterogeneous . The BVP 
satisfies the maximum principle and its solution stays between and 2. 

_1 (see Fig. 



Uniform meshes in metric 



(see Fig. 



10 



10 



'b)) are used 



) ^a)) and in metric 6 (eh) 

for this example, where 6 (eh) is a scalar function depending on the error estimation e^. It is known 
|19| that for these meshes, the linear finite element solution of the anisotropic diffusion problem 
without convection and reaction terms satisfies DMP. 



Contours of linear finite element solutions are shown in Figs. 11 and 12 One can see that for 
both types of meshes, undershoots occur for a relatively coarse mesh and vanish for a finer mesh. 
The result is consistent with the theoretical prediction in the previous section. □ 



6 Conclusions 



In the previous sections we have developed a mesh condition (25) under which the linear finite 
element solution defined in ^ for the general anisotropic diffusion problem ([!]) and ^ involving 
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(a): N = 2312 



(b): N = 4232 





Figure 8: Contours of linear finite element solutions for Example 5.3 
(a): b = [200, 200] T (b): N = 3200 



N: numbers of elements 



Figure 9: The undershoot, —u m i n , is plotted as a function of the number of elements N in (a) and 
as a function of 



in (b) for Example 5.3 



convection and reaction terms to satisfy a discrete maximum principle. Loosely speaking, the 
condition requires that the dihedral angles of all elements of the mesh be Odl&Hoo/i+Hclloo/i 2 )— acute 
when they are measured in the metric specific by the inverse of the coefficient matrix, where b and 
c are the coefficients of the convection and reaction terms, respectively. Moreover, we have shown 

an O 



that in two dimensions a weaker condition, (34) or (36) 



folloo/i + || c|| oo/i 2 ) perturbation 
of the generalized Delaunay condition developed in [12J, is sufficient for the linear finite element 
solution to satisfy a discrete maximum principle. Finally, it is worth pointing out that many existing 
mesh conditions such as those developed in Ciarlet and Raviart [6] (for isotropic diffusion without 
convection terms), Strang and Fix [27] (for 2D isotropic diffusion without convection terms), Wang 
and Zhang [28] (for isotropic diffusion with convection and reaction terms), Li and Huang [19] (for 
anisotropic diffusion without convection and reaction terms), and Huang [12] (for 2D anisotropic 



diffusion without convection and reaction terms) are special cases of mesh condition ( 25 ) or ( 34 ) 
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(a): N = 2050 



(b): N = 1822 





Figure 10: Uniform meshes in metric D 1 (a) and in metric 0(e^)D 1 (b) used for Example 5.4 



fa): N = 2050 



(b): N = 4184 





Figure 11: Contours of linear finite element solutions obtained with uniform meshes in metric 



for Example 5.4 



10931004 and 40906048. 
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